home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATH1 / NLIN3.PAS < prev    next >
Pascal/Delphi Source File  |  1985-04-03  |  4KB  |  193 lines
  1. { die ln-funktion bei MT+ ist fuer Werte kleiner 1E-5 aeussertst langsam,
  2.   so dass man das Gefuehl bekommt, das dass Programm sich aufhaengt. In-
  3.   folgedessen konnte ich dieses Programm nicht verifizieren. -Juergen }
  4.  
  5. program nlin3;        { -> 310 }
  6. { Pascal program to perform a nonlinear least-squares fit for the diffusion
  7.     of Zn in CU }
  8.  
  9. const    maxr    = 20;        { data prints }
  10.     maxc    = 4;        { polynomial terms }
  11.     r    = 1.987;    { gas constant }
  12. type
  13.     index    = 1..maxr;
  14.     ary    = array[index] of real;
  15.     arys    = array[1..maxc] of real;
  16.     ary2    = array[1..maxr,1..maxc] of real;
  17.  
  18. var
  19.     x,y,y_calc    : ary;
  20.     t,d,ex        : ary;
  21.     coef        : arys;
  22.     i,n        : integer;
  23.     nrow,ncol    : integer;
  24.     done,error    : boolean;
  25.     correl_coef,srs,
  26.     a,b,x2        : real;
  27.  
  28. external procedure cls;
  29.  
  30. procedure get_data(var x,y: ary;
  31.            var n: integer);
  32. { get values for n and arrays t,d }
  33.  
  34. var    i    : integer;
  35.  
  36. begin
  37.   n:=7;
  38.   t[1]:=600.0;    d[1]:=1.4E-12;
  39.   t[2]:=650.0;    d[2]:=5.5E-12;
  40.   t[3]:=700.0;    d[3]:=1.8E-11;
  41.   t[4]:=750.0;    d[4]:=6.1E-11;
  42.   t[5]:=800.0;    d[5]:=1.6E-10;
  43.   t[6]:=850.0;    d[6]:=4.4E-10;
  44.   t[7]:=900.0;    d[7]:=1.2E-9;
  45.   for i:=1 to n do
  46.     begin
  47.       x[i]:=1.0/(t[i]+273.0);
  48.       y[i]:=d[i]
  49.     end
  50. end;        { proceddure get data }
  51.  
  52. procedure write_data;
  53. { print out the answers }
  54. var    i    : integer;
  55. begin
  56.   writeln;
  57.   writeln;
  58.   writeln('     I      TC     D      DCALC');
  59.   for i:=1 to n do
  60.     writeln(i:3,t[i]:8:0,d[i],' ',y_calc[i]);
  61.   writeln; writeln(' Coefficients ');
  62.   writeln(coef[1],' constant term');
  63.   for i:=2 to ncol do
  64.     writeln(coef[i]);        { other terms }
  65.   writeln;
  66.   writeln('D0=',a:7:2,' cm sq/sec.');
  67.   writeln('Q =',(-r*b/1000.0):8:2,'kcal/mole');
  68.   writeln;writeln('SRS= ',srs:8:4)
  69. end;        { write_data }
  70.  
  71. procedure func(b: real;
  72.       var fb,dfb: real);
  73. var    i        : integer;
  74.     s1,s2,s3,s4,s5,s6,
  75.     ex1,ex2,xi,
  76.     x2,yi,y2    : real;
  77. begin
  78.   s1:=0.0;
  79.   s2:=0.0;
  80.   s3:=0.0;
  81.   s4:=0.0;
  82.   s5:=0.0;
  83.   s6:=0.0;
  84.     for i:=1 to n do
  85.       begin
  86.     xi:=x[i];
  87.     x2:=xi*xi;
  88.     yi:=y[i];
  89.     y2:=yi*yi;
  90.     ex1:=exp(b*xi);
  91.     ex[i]:=ex1;
  92.     ex2:=ex1*ex1;
  93.     s1:=s1+xi*ex2/y2;
  94.     s2:=s2+ex1/yi;
  95.     s3:=s3+xi*ex1/yi;
  96.     s4:=s4+ex2/y2;
  97.     s5:=s5+2.0*x2*ex2/y2;
  98.     s6:=s6+x2*ex1/yi
  99.       end;
  100.     fb:=s1*s2-s3*s4;
  101.     dfb:=s2*s5-s1*s3-s4*s6;
  102.     a:=s2/s4
  103. end;        { func }
  104.  
  105.  
  106. procedure newton(var x: real);
  107. const    tol        = 1.0E-6;
  108.     max        = 20;
  109. var    fx,dfx,dx,x1    : real;
  110.     i        : integer;
  111.  
  112. begin    { newton }
  113.   error:=false;
  114.   i:=0;
  115.   repeat
  116.     i:=i+1;
  117.     x1:=x;
  118.     func(x,fx,dfx);
  119.     if dfx=0.0 then
  120.       begin
  121.     error:=true;
  122.     x:=1.0;
  123.     writeln('ERROR: slope zero')
  124.       end
  125.     else
  126.       begin
  127.     dx:=fx/dfx;
  128.     x:=x1-dx;
  129.       end
  130.   until
  131.     error or
  132.       (i>max) or
  133.     (abs(dx)<=abs(tol*x));
  134.   if i>max then
  135.     begin
  136.       writeln(chr(7),'ERROR: no convergence in ',max,' loops');
  137.       error:=true
  138.     end
  139. end;        { newton }
  140.  
  141. procedure nlin(x,y:    ary;
  142.     var y_calc:    ary;
  143.            n:    integer);
  144. { fits the diffusion equation through n sets of x and y pairs of points }
  145. var
  146.     resid        : ary;
  147.     matr        : ary2;
  148.     i        : integer;
  149.     xi,yi,sum_x,
  150.     sum_y,sum_y2,b1,
  151.     sum_xy,sum_x2    : real;
  152. begin        { nlin }
  153.   ncol:=2;    { number of terms }
  154.   sum_x:=0.0;
  155.   sum_y:=0.0;
  156.   sum_xy:=0.0;
  157.   sum_x2:=0.0;
  158.   for i:=1 to n do
  159.     begin
  160.       xi:=x[i];
  161.       yi:=ln(y[i]);
  162.       sum_x:=sum_x+xi;
  163.       sum_y:=sum_y+yi;
  164.       sum_y2:=sum_y2+yi*yi;
  165.       sum_xy:=sum_xy+xi*yi;
  166.       sum_x2:=sum_x2+xi*xi
  167.     end;
  168.   b:=(sum_xy-sum_x*sum_y/n)/(sum_x2-sqr(sum_x)/n);
  169.   newton(b);
  170.   coef[1]:=a;
  171.   coef[2]:=b;
  172.   srs:=0.0;
  173.   for i:=1 to n do
  174.     begin
  175.       y_calc[i]:=a*ex[i];
  176.       if y[i]<>0.0 then
  177.     resid[i]:=y_calc[i]/y[i]-1.0
  178.       else resid[i]:=y[i]/y_calc[i]-1.0;
  179.       srs:=srs+sqr(resid[i])
  180.     end
  181. end;         { nlin }
  182.  
  183.  
  184. begin        { main program }
  185.   cls;
  186. writeln(' start get_data ');
  187.   get_data(x,y,n);
  188. writeln(' start nlin ');
  189.   nlin(x,y,y_calc,n);
  190. writeln(' start write_data ');
  191.   write_data
  192. end.
  193.